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Abstract 

In the field of neutrino astronomy, large volumes of optically transparent mat- 
ter like glacial ice, lake water, or deep ocean water are used as detector media. 
Elementary particle interactions are studied using in situ detectors recording time 
distributions and fluxes of the faint photon fields of Cherenkov radiation generated 
by ultra-relativistic charged particles, typically muons or electrons. 

The Photonics software package was developed to determine photon flux and 
time distributions throughout a volume containing a light source through Monte 
Carlo simulation. Photons are propagated and time distributions are recorded through- 
out a cellular grid constituting the simulation volume, and Mie scattering and ab- 
sorption are realised using wavelength and position dependent parameterisations. 
The photon tracking results are stored in binary tables for transparent access 
through ansi-C and C++ interfaces. For higher-level physics applications, like simu- 
lation or reconstruction of particle events, it is then possible to quickly acquire the 
light yield and time distributions for a pre-specified set of light source and detector 
properties and geometries without real-time photon propagation. 

In this paper the Photonics light propagation routines and methodology are 
presented and applied to the IceCube and Antares neutrino telescopes. The way 
in which inhomogeneities of the Antarctic glacial ice distort the signatures of ele- 
mentary particle interactions, and how Photonics can be used to account for these 
effects, is described. 
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1 Introduction 



In optical high energy neutrino astronomy light from particle physics events 
is observed using a large number of detectors placed deep in glacial ice or in 
ocean or lake water. Successful simulation and reconstruction of such events 
relies on accurate knowledge of light propagation within the detector medium. 
Light propagating through even the clearest water or ice is affected by scat- 
tering and absorption. For light sources and receivers separated by distances 
comparable to the photon mean free path, scattering effects can neither be 
analytically calculated nor ignored. The typical scattering lengths in these de- 
tection media are tens to hundreds of metres. Since this scale is comparable 
to the typical detector separation, detailed simulation of the photon propa- 
gation is required to obtain information necessary for event simulation and 
reconstruction. The problem is complicated further by the anisotropy of the 
light emitted in particle interactions and the heterogeneity of detector media. 
Photonics is a freely available software package [1] containing routines for 
detailed photon Monte Carlo simulations, which take into account such com- 
plexities to provide, in tabulated form, the photon flux distribution throughout 
a specified medium for an input light source. 

With Photonics the photon flux and time distributions can be tabulated for 
an arbitrarily large volume of a propagation medium, for a user defined range 
of light source and detector properties. This means that once a Photonics 
table set has been generated for a class of light sources and detectors, it is pos- 
sible to quickly and transparently acquire the light yield and time distributions 
without any need for real-time photon propagation during, for example, par- 
ticle physics event simulation or reconstruction. This is made possible by the 
Photonics table reader library, with which a user (program) can dynamically 
query the pre-calculated tables by specifying the locations and geometrical re- 
lations between light sources and detectors. A simulation chain for a complete 
experimental setup can be achieved by using these interfaces and applying 
detector specific details such as modelling of electronics, data acquisition, and 
triggers. For event reconstruction, PHOTONICS provides probability density 
functions for arrival times of independent photons and the expected number 
of detected photons. 

In this paper we first introduce the relevant physics of the photon propagation 
(Section 2) and the details of the Photonics implementation (Section 3). We 
then compare our results with observations of calibration light sources in sea 
water and glacial ice (Section 4). In Section 5, we present some photon track- 
ing results relevant to the detection of neutrinos with the IceCube neutrino 
telescope. 
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2 Light propagation in diffuse media 

Our goal was to model the transport of light through glacial ice and water. 
Photon propagation depends on the optical properties of the medium, in par- 
ticular on the velocity of light and the absorption and scattering cross sections. 
Glacial ice is optically inhomogeneous because of depth dependent variations 
in temperature, pressure, and concentrations of air bubbles and insoluble dust. 
Since the dust deposits track climatological changes and are therefore assumed 
to be arranged in horizontal layers, their effect is parameterised as a vertical 
variation of the optical properties. In addition to this spatial variation, the 
wavelength dependence of the medium parameters must be taken into ac- 
count. Before describing the implementation to achieve our goal, we review 
the optical quantities that must be considered in the simulation, using nota- 
tion in which the wavelength dependence is left implicit. 

The time of light travel is determined by the group velocity of light, which is 
given by the group refractive index n g , while various transmission and scat- 
tering coefficients depend on the phase velocity [2] and its index of refraction 
n p . 

Absorption of visible and near UV photons in pure water and ice is due to 
electronic and molecular excitation processes and is characterised by the ab- 
sorption length A a . Measurements of light attenuation have been performed 
at relevant wavelengths in lake water by the Baikal [3] collaboration, and in 
sea water by the Dumand [4], Nestor [5], Antares [6], and Nemo [7] col- 
laborations. The Amanda collaboration has developed an empirical model 
for optical absorption in deep glacial ice by combining laboratory and in situ 
measurements [8]. 

Photon scattering by scattering centres of general sizes is described by Mie 
scattering theory [9], which for any wavelength and scattering centre size 
gives the scattering angle distribution, the phase function. Rahman scattering, 
where the scattering centre is affected by the scattered photon, and Brillouin 
scattering, where photons are scattered on (thermal) density fluctuations, may 
result in a change of photon energy. However, these processes are subdomi- 
nant to Mie scattering in both glacial ice [10] and sea water [6] where light 
is scattered by centres of very different types and sizes: from ice crystal point 
defects to air bubbles and mineral grains in ice, and from biological matter to 
sediment particles in water. 

In natural ice and water it is difficult to determine the phase function from in 
situ measurements. Instead, the Amanda and Antares collaborations have 
used calibration light sources to determine the propagation characteristics 
assuming certain forms of the scattering angle distributions. In the case of 
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ice, a one parameter Henyey-Greenstein (HG) phase function is often used, 
approximating Mie scattering under the assumption that scattering is forward 
peaked [11]. For water, a two parameter phase function is more useful. For this 
paper we mostly use the single parameter HG phase function 

1 -r 2 

PHG(cOsfl) = — — - - ¥ , (1) 

2(1 + T 2 — 2r COSt^ja 

which is completely characterised by the r parameter, the mean of the cosine 
of the scattering angle 0, 

t = (cos 6) = / £>hg(cos#) cos#d(cos#). (2) 



The absorption length, A a , and scattering length, A s , are the mean free paths of 
exponential distributions. The probability density function for the path length 
s to the next scatter is 

f M - dF(s) e ~ S/As (*\ 



where F(s) is the probability distribution function. 

When determining ice or water optical properties, there is a degeneracy be- 
tween A s and r. One therefore considers the effective scattering length, A c , 
defined as 



which in anisotropic scattering is analogous to the (geometric) scattering 
length A s in isotropic scattering; it is the distance which light propagates 
through a turbid medium before the photon directions are completely ran- 
domised. Consider a collimated light pulse injected into a non-absorbing medium. 
In this case, the photons are on average scattered at successive steps of length 
A s and the projection of the net velocity vector on the original direction is 
decreased on average by r = (cos 6) in each scattering step (in which all pho- 
tons are scattered) [12]. Hence the injected light is effectively transported a 
forward distance of 

A s ]>>^-^ = A e , (5) 



and A c has a natural interpretation as the distance that the centre of gravity 
of the photon cloud advances, in the limit of many scatters. 
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3 Monte Carlo simulation implementation 

In this section we present the main ingredients in our photon Monte Carlo 
simulation implementation. The end product is a set of photon flux density 
tables describing the evolution of the light field around a user-defined source. 
For a given light source, a user-specified large number of photons is generated 
according to the source characteristics. The photons are then tracked and their 
contribution to the overall light field is determined and recorded in a cellular 
grid throughout a user defined portion of space, the simulation volume. The 
sensor locations are not fixed, but are dynamically specified when accessing 
the simulation results. The photon intensity and time distributions are stored 
in a six dimensional binary table. Four of these dimensions are for the spatial 
and temporal location in the simulation volume with respect to the emission 
point. As the acceptance of the light sensors is assumed to be azimuthally 
symmetric, around the vertical axis in heterogeneous media and around any 
axis in homogeneous media, the photon impact direction is characterised by 
the zenith angle alone, constituting the fifth dimension. The sixth dimension 
is the angle from the light source principal axis at which a photon is emitted 
(this dimension can be used, for example, to reweight the flux tables for a 
different emission profile). These latter two dimensions are usually integrated 
over when the photon tables are used with detector simulations. In this case the 
wavelength and angular sensitivity of the detector elements need to be folded 
in as the photons are recorded, using recording weights which can be specified 
in functional or tabular form. For each light source position and orientation 
one table is produced. A set of tables describes a range of source locations 
and directions, valid for the specified class of light sources and detectors. 



3. 1 Media and light source properties 

The parameters describing the optical medium (n g , n p , r, A e , A a ) are taken to 
be functions of wavelength and a spatial dimension Z, typically specifying the 
ice or water depth. Thus, the propagation medium is divided into horizontal 
regions which differ by their optical properties. For media where the single 
parameter HG approximation does not provide an adequate description of the 
scattering it can be replaced by other phase functions. An example of this is 
the treatment of sea water in Section 4.1. 

A single simulation run begins with the injection of a photon with wavelength 
and emission direction chosen from user-specified probability distributions and 
at a user-specified location. Our procedures support point-like (infinitesimal) 
light sources, where all emitted photons originate from the same point, and 
volume light sources where the emission is distributed over a volume. An ex- 
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ample light source, essential for neutrino astrophysics, is that of a Cherenkov 
emitter. In PHOTONICS a Cherenkov emitter is a point-like light source with 
a Cherenkov wavelength spectrum and an angular emission in a Cherenkov 
cone [13]. In this case the emission is azimuthally symmetric around the prin- 
cipal axis of the Cherenkov emitter. Closely related are light sources composed 
of many short Cherenkov emitting tracks, such as electromagnetic cascades. 
Another category of point-like sources are laser or LED light sources, and in 
Section 4 we compare our simulation of such sources with observations. Con- 
tinuous and extended light sources are composed by integration of infinitesimal 
sources. For example, the light distribution due to a relativistic muon is com- 
posed by integrating a series of infinitesimal Cherenkov emitters over space 
and time, as described in Section 3.4. Simulation results for both point-like 
and extended Cherenkov emitters are presented in Section 5. 



3.2 Coordinate systems for photon flux recording 

The Cherenkov light source example possesses cylindrical symmetry. Although 
this symmetry is typically broken by the response of the propagation medium 
and the detector, a cylindrical coordinate system aligned with the source's 
symmetry axis is a natural choice for the recording of the light flux, and it 
is therefore used in the following. In addition to cylindrical (p, I, <p) coordi- 
nates, we have also included functionality to allow the flux to be recorded in 
spherical (r, 0) or Cartesian (e 1 ,e 2 ,e 3 ) coordinates. A grid over the spatial 
coordinates defines cells in which the photon flux density and time distribu- 
tions are averaged. The user-specified region of space covered by these cells 
constitutes the recording volume. 

In some dimensions it can be desirable to have denser binning close to the 
origin. Therefore, uniform as well as linearly increasing bin sizes are supported 
for p, I, and the time t. In addition to specifying the spatial coordinates where 
the light flux is recorded relative to the source, the location of the source and 
the direction of its axis of symmetry with respect to the medium are needed. 
These can be characterised using just two coordinates, the depth Z s and the 
zenith angle S . This is because of the horizontal symmetry of the medium, 
as discussed in Section 3.1. Fig. 1 shows the coordinates and a recording cell 
in which the average flux is recorded. 

If in addition to the medium symmetry around z we assume that the light 
source is symmetric with respect to the sign of <j) (besides any light source 
that is azimuthally symmetric around its axis this is the case for some non- 
isotropic LED calibration light sources used in glacial ice [8,14]), the flux 
needs only to be tabulated for between 0° and 180°. For heterogeneous 
media we require that the sensor acceptance and the medium properties are 
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symmetric around the same axis z. This is not required for homogeneous 
media since the coordinate system can be rotated make the sensor symmetry 
axis collinear with the z axis used to define source location and orientation. 
Photonics can be extended to handle cases where the two symmetry axes do 
not align in heterogeneous media, by adding further dimensions to the photon 
tables. 



Dimension 




Bins 


Low 


High 


Radial 


P 


30 


m 


500 m 


Longitudinal 


I 


51 


-500 m 


500 m 


Azimuthal 


4> 


10 


0° 


180° 


Time 


t 


50 


ns 


6000 ns 



Table 1 

An example of table binning using cylindrical coordinates. The rotational symmetry 
of the emitter and the horizontal symmetry of the medium implies an azimuthal 
symmetry in <f>, so that the flux is the same at —(f) and <fi. 

A binning example is shown in Table 1 for a cylindrical coordinate system. As 
each single precision floating point number requires four bytes, the size of this 
example table would be 3 megabytes. To get the total table set size this must 
be multiplied with the number of light source positions Z s and zenith angles 
6 S of interest. For example, with 50 source depths and 20 source angles, the 
total table set size is 3 gigabytes. 

3.3 Photon propagation 

Starting at their point of origin, the photons are tracked and recorded along 
straight-line paths between successive scattering points using one of two meth- 
ods. In the volume-density mode, photons are recorded at equidistant record- 
ing points along their paths. In the area-crossing mode, photons are recorded 
at every surface-crossing into a cell. In the latter mode, the propagation step 
is dynamically calculated to bring the photon to its next cell boundary. The 
photon's propagation properties are always updated at medium boundaries 
and scattering points, regardless of recording method. 

When a scattering point is reached, the photon's direction is changed by an 
angle randomly drawn from the selected phase function, typically the HG 
function in Eq. (1), and in a uniform random azimuthal direction. At this 
point, the distance to the next scattering point is drawn from an exponentially 
distributed random variable with a mean value of the local scattering length. 

Photon absorption is taken into account by successively updating the photon 
survival weight w as the photons propagate through regions of different ab- 
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Figure 1. The recording cell geometry with variables used for binning of photon 
flux data. Photons are emitted from a user-defined light source at Z s and their flux 
and time distributions are recorded and averaged in all spatial cells the photons 
traverse (one of which is shown as a shaded volume). These cells are defined in 
a coordinate system aligned with the source's principle axis I, which is tilted by 
7r — S with respect to the medium symmetry axis z. The angle 4> is defined to be 
zero where the radial vector is maximally aligned with z. The azimuthal direction 
<3? s is degenerate since the medium properties are assumed to be symmetric around 
the z axis. 

sorption. For a photon tracked n steps, each within a locally homogeneous 
medium, the weight is given by 



where As« is the length of step i in a region with absorption length X a ^. The 
weight is updated every time the photon is scattered, enters a new medium 
region, or is recorded. 

When a photon enters a medium region with different scattering and ab- 
sorption parameters, these are updated, from (A s , A a ) to (\' s ,X' a ), at the re- 
gion boundary. At this point, the remaining distance to the next scatter is 
s' = s\' s /X s , where s would have been the remaining distance to the next scat- 
ter in the former region with scattering length A s . Refraction at the boundary 
is supported but reflection is ignored since it is assumed that refractive index 
variations are continuous. 

During its propagation, the flux contributed by a photon is recorded either in 




(6) 
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each spatial cell it enters (in area-cossing mode) or each time it completes a 
propagation step (in volume-density mode). The photon flux <3> (particles per 
area and time) at any point (p, I, <fi) at a time T after the emission from a light 
source at a depth Z s pointing in a direction O s is denoted &(Z S , 6 S , p, l,(f>,T — 
to(p,l)), where to is a reference time typically chosen to be the first time 
causally connected to the light emitted by the source, 

i o (p,0 = — V^+^ (7) 



where c is the speed of light in vacuum and n rc f is a user-specified reference 
refractive index. This reference time convention is appropriate for point-like 
stationary light sources only. For fast moving light sources such as muons a 
different expression is used which takes into account the more complicated 
causality condition. The residual time t = T — t is the time delay caused by 
scattering, relative to the propagation time for a photon travelling in a straight 
line. Photons are tracked until their residual times exceed a user-specified 
value. The tracking of a photon is terminated if it leaves the simulation volume 
(which can be larger than the recording volume to allow the photon to scatter 
back into the recording volume) or if its survival weight drops below a pre-set 
value. 

The probability density function for a photon flux at time t is given by 

f (7 ft n 1 rh A - $(^s,e s ,p,Z,<M) 

j P df (A, Ws, p, t, (p, t) - ttt~t\ — mv - ' y°> 

where I is the time integrated photon flux, or intensity, 

oo 

J(Z s ,6 s ,p,/,0) = J $(Z B ,G B ,p,l,<l>,t)dt. (9) 

— oo 



Since photon fluxes are additive we can determine the time distribution of a 
combination of light sources through 

Mt) = EdMl = M (10) 



The way in which the photon intensity in a spatial cell is calculated depends on 
the recording method. In the area-crossing method, the contribution of each 
photon to the total flux in the cell depends on the projected surface area of 
the cell as seen in the direction of the contributing photon as it crosses the cell 
boundary. In the volume-density method, photons contribute at equidistant 
recording points along their paths, so that the contribution is proportional 
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to the number of recording points that fall in a given cell. The respective 
equations for the calculation of the observed intensity, per emitted photon, 
are 



AreaCrossing 



1 N 



(11) 



7=1 



1 N «7 



As 



VolumeDensity 



• v 7 =lfc=l 



(12) 



where the 7 sum runs over the N simulated photons, A± is the recording cell 
area projected perpendicularly to the direction of the photon, V is the volume 
of the cell, and As is the recording point separation along the photon path in 
the volume-density method. The quantity n 7 is the number of recording points 
inside the cell for a given photon. Hence, in the area-crossing mode, photons 
are recorded at every surface-crossing into a cell, whereas in the volume-density 
mode, they are recorded at equidistant recording points. The photon weight u> 7 
is a product of the absorption-induced survival weight, Eq. (6), and optional 
user-specified sensor detection efficiencies. 

The calculation and recording of a photon's contribution to the flux in a 
recording cell is computationally expensive, but a suitable choice of method 
(area-crossing or volume-density) can speed up the simulation. To optimise for 
speed one compares the step size As with the scale dimension of a recording 
cell -Dccii- If As < -Dceii, the area-crossing method is competitive; otherwise the 
volume-density method should be used. Hence the area-crossing method can 
result in faster code execution for large detection volumes with large recording 
cells, while the volume-density method is faster for small, dense recording cell 
configurations. To ensure unbiased sampling in volume-density mode even for 
As 3> -Dceii, the path length to the first recording point after emission is drawn 
from a uniform distribution between and As. To maximise the number of 
independent sampling points for a given execution time, As is balanced against 
the total number of photons N. A large As allows a large N, but if too large 
the number of recordings per execution time will drop as more time is spent 
propagating photons between recording points. Using As A s is often a good 
compromise. 

Another way to improve the speed of the Monte Carlo simulation is through 
importance-weighted scattering. This means that the photons are propagated 
using scattering parameters which are different from those of the scattering 
situation at hand in order to get higher statistics at low probability phase 
space locations. When solving a Monte Carlo problem for random numbers xi 
(for example, scattering angles) from a probability density distribution f\(x), 
we can choose to instead sample from another distribution /^(a?) while apply- 
ing a weight of fi(xi) / fzixi)- As an example, a user might want to oversample 
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straighter paths to enhance the statistics for early photons at large source- 
receiver distances. The user could then simply scale the scattering length by 
a factor k > 1, propagate photons with scattering described by fk\ s ( s ) ( see 
Eq. (3)), and reweight each photon by f\ 3 (s)/fk\ 3 (s) for every scatter it ex- 
periences. Alternatively, the mean cosine of the scattering angle, r, can be 
modified to achieve a similar effect using Eq. (4). The effective scattering 
length is made a factor k longer by r' = (k — 1 + r)/k. The corresponding 
weight is fM/fAO). 

We have described how the photon flux density is calculated in cells throughout 
the simulation volume. In neutrino astronomy, the quantity of interest is the 
number of photons detected by an optical sensor such as a photomultiplier. 
Since the sensor response is usually wavelength and angle dependent but we 
do not want to store the wavelengths and arrival directions of the simulated 
photons, Photonics includes the option to fold the sensor response into the 
simulation when the tables are generated. When this option is selected, user- 
supplied wavelength and angular efficiency files are used to weight the photon 
flux density to obtain the detected photon flux. 



3.4 Propagating light sources 



To this point we have discussed how to obtain a set of tables describing the 
photon fluxes for a range of locations and orientations of a point-like, station- 
ary light source. A propagating light source can not be satisfactorily approxi- 
mated as a flash of light from a single spatial point. We discuss in the following 
the modelling of Cherenkov light from high energy muons, which give rise to 
kilometre scale Cherenkov emitting tracks in water or ice. Our method, how- 
ever, applies also to other line-like light sources such as high energy tauons. 

To calculate the photon flux distribution generated by a muon we integrate 
over the photon flux distributions of many point-like Cherenkov emitters with 
O s given by the muon direction. A set of point-like photon tables provide the 
differential light flux $ po int at any space-time location from sources at any 
causally connected location. This serves as integration kernel, 

(-*point(>£s(^s)> £C rec ,t rec ) = 7T - ^point (•Es(^s) i -^rec ; ^rec) ; (13) 



where (x s ,t s ) and (x rec ,t Tec ) are the emission and receiver coordinates. The 
light distribution of a propagating muon is thus generated by convolving this 
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kernel with the track of the muon, x^t^), so that 

*Mstop 



Photonics provides the capability to efficiently perform this integration for 
a series of fixed muon directions and locations and to store the resulting light 
fluxes in tables like those for the point-like light sources. These tables are then 
used to obtain the light flux for muon tracks through any part of the simulation 
volume by interpolation (Section 3.5). Although the construction of muon 
events through the integration in Eq. (14) can in principle also be done event 
by event in a detector simulation or event reconstruction, such an approach 
would often be significantly slower since the number of considered events, and 
therefore the number of required flux calculations, is typically larger than the 
number of fixed muon directions and locations needed to adequately cover the 
relevant range. 

To provide the photon flux for any given track length without having to dy- 
namically perform the time consuming integration of Eq. (14), we have de- 
veloped a scheme based on the subtraction of semi-infinite tracks. The flux 
at (aj rec , tree) arising from a finite muon starting at Xa and stopping at x-q 
can be expressed as the difference between two semi-infinite tracks. For two 
semi-infinite starting tracks, one (denoted /xa->) starting at the point x& and 
the other (//b-0 starting at x B located further along the muon track, we write 

*^Vab (*^rec; tree) ^*/Ua— ► (^rec? ^rcc) ^Vb-> (^rec; tree) j (^^) 



and the probability density function can be written 

J MAB ~ t _ t \ lV J 



using Eq. (8). This construction of finite tracks is depicted in Fig. 2. Consider 
the point p , close to the finite track start point. At this point, the contribution 
from $ MB ^ is comparably small, so that $ MAB ~ ^a^- However, for the point 
p l close to the end of the finite track, $ MA ^ ~ $ MB ^ and the calculation of 
^Vab becomes sensitive to exact numerical cancellation. It is then numerically 
superior to subtract two semi-infinite stopping tracks, /x_»a and //_»b, and write 
^Vab = ^V^b — ^V^a- O ur algorithm dynamically chooses between these two 
descriptions depending on which of the endpoints is closer to the observation 
point. 

For applications in which tracks can be regarded as completely infinite, the 
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Po 



II 



III 



IV 



Figure 2. Schematic view of flux tables for four collinear tracks (offset in x for 
clarity). The light distribution of a finite track (I) is the difference between two 
semi- infinite tracks (II) and (III). The dotted rectangles show the outer limits of 
the semi-infinite tables. The gray lines represent example isointensity contours. The 
representation of photon flux at a point p 2 far from the starting and stopping points, 
can be done in two ways. Either, as in (III), by considering the edge of the best 
matching semi-infinite table, or, as in (IV), with an infinite track table with only a 
single I bin. 



flux tables are made smaller by removing redundant information. At a given 
observation point in the medium, the flux from an infinite muon track with 
table origin Zi, giving the point a longitudinal coordinate l\, is identical to the 
flux in a table with origin Z 2 where the point is at h = h + (Z2 — Z\)j cos O s . 
Therefore in each infinite muon table only one I bin is retained, typically at 
1 = 0. Fig. 2 illustrates both ways of accessing information about virtually 
infinite tracks: with infinite tables and with semi-infinite tables. 



The size of a set of infinite tables is typically about 1% of the size of the corre- 
sponding semi-infinite description, which can be tens of gigabytes. Hence the 
infinite tables can easily be loaded into computer memory all at once, which 
is particularly useful in event reconstruction where it can be hard to estimate 
in advance the properties of an event to be reconstructed. Reconstruction 
and large table support is discussed in the following section, describing the 
Photonics reader library. 
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3.5 Using the photon flux tables for event simulation and reconstruction 

For event simulation and reconstruction, the photon flux tables are accessed 
using either a set of ansi-C procedures, or a more abstract (root compliant) 
user interface written in C++, both provided with the Photonics package. 

The full simulation of for example an ultra relativistic muon crossing the 
detector volume is performed by first propagating the muon through the de- 
tector medium with a charged particle propagator such as MMC [15]. This 
results in a list of light generating subevents such as minimum ionising muon 
track segments and associated electromagnetic showers induced by stochas- 
tic energy losses. The detector specific simulation program can then query 
the corresponding Photonics table information to obtain the number and 
time distribution of detected photons in each detector module from each such 
subevent. Photonics comes with a variety of idealised light emission profiles, 
such as those of minimum ionising muons and point-like electromagnetic and 
hadronic showers. 

A set of tables is loaded into memory using the table reader library. The 
user can then query the photon flux tables by specifying the location (x s ) 
and orientation (G s , <3> s ) of the source, and the location of the light detectors 
(iCrec). Two additional source characteristics are optional in the query: the 
source length L (applicable to finite muons) and the energy E used to scale 
the light source intensity. In an experiment simulation, the user first requests 
the expected number of detected photons N(x s , S , $ s , x TCC , E, L). This query 
also returns a table reference which can be used to get photon arrival times 
randomly drawn from the tabulated time distribution at the corresponding 
coordinates. For event reconstruction, Photonics also provides the arrival 
time probability density function / pd f (x s , S , <J> S , x Tec , L, t), which can be used 
by track-fitting algorithms, for example maximum-likelihood routines. Both 
/pdf and the photon intensity (giving N) are naturally continuous in L and 
E, and are made continuous in x s , 6 S , $ s , x Tec , and t by multidimensional 
linear interpolation. The interpolation of time distributions is flux weighted, 
in agreement with Eq. (10), 

/pdf = ^' UtIl j\ with J2i^i = 1, (17) 

where u)i are the interpolation weights. When the flux for a requested source 
direction @ s is interpolated from two surrounding tabulated directions 0i and 
G>2, these tables are first (implicitly) rotated to @ s . The receiver coordinates 
then identical for the 0i and O2 table. An analogous approach is used 
for the source origin x s . 
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It is sometimes necessary to convolve the photon time distributions with the 
detector time response function or the emission time profile. This can be done 
with a provided routine operating on the photon flux tables. Convolution 
with a Gaussian or with one of two light source time distributions with longer 
positive tails [6] (see Section 4) has been implemented. 

Any number of table sets can be loaded simultaneously (limited only by avail- 
able memory), making it possible to simulate or reconstruct the cumulative 
signal from different types of light sources (and detectors). This is useful for 
example when simulating muons with secondary electromagnetic showers from 
a primary muon track or several coincident muons. Since detailed photon ta- 
ble sets are often many times larger than the primary memory of a computer, 
we provide some memory management tools. Users can dynamically load and 
unload tables, and select loading of tables corresponding to limited ranges of 
depths (Z s ) and light source angles (6 S ). In experiment simulations this is 
particularly useful since the parameters for the simulated particles are known. 
In event reconstruction, the loaded tables should cover the phase space taken 
into account in the fitting algorithms. If memory limitations restrict this phase 
space, reconstruction can be performed for subregions, on events predeter- 
mined to lie within subregions small enough for the corresponding photon 
tables to fit in memory. Such preselection can be done with a set of more 
coarsely binned Photonics tables or with other first-guess approaches (for 
muons, the much smaller tables of the infinite description can be used, see 
Section 3.4). 

In addition to what has already been mentioned, the Photonics package in- 
cludes several other tools for processing of the photon flux tables, including 
integration in any dimension and conversion between differential and cumula- 
tive time distributions. 



4 Comparison with observations 

In this section we use measurements with artificial calibration light sources 
in sea water and glacial ice to demonstrate that Photonics reproduces the 
general behaviour of the observed photon distributions. 

4-1 Modelling of natural water and Antares ' Mediterranean water surveys 

The Antares collaboration is constructing a 0.1 km 2 water neutrino detector 
in the deep Mediterranean sea [6]. In the present design the detector has 12 
vertical strings, each of which has an instrumented height of about 350 m 
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and consists of 25 storeys with three optical sensors each. Three strings were 
deployed in 2006 and the remaining strings are scheduled to be installed during 
2007. 



Absorption and effective scattering lengths in the deep Mediterranean sea wa- 
ter have been investigated by the Antares collaboration [6]. The surveys 
were performed during several seasons using a calibrated setup of isotropic 
light sources, one in the blue at 473 nm and one in the UV at 375 nm. For 
blue (UV) a A a of 60 (26) m and a A e of 265 (122) m with 15% time vari- 
ability are quoted by Antares. The details of the experimental setup, such 
as the light emission time profile of the source and detector efficiencies, play 
a large role for the photon flux time distribution profile in media, like ocean 
water, where scattering is weaker than absorption. In addition, since the light 
is typically observed at distances shorter than or comparable to the effec- 
tive scattering length, the scattering phase function plays an even larger role. 
Antares assumed a weighted sum of a Petzold and a molecular (Einstein- 
Smoluchowski) distribution. The Petzold distribution has r = 0.924, and is 
here approximated by a HG distribution, while the molecular distribution is 
approximated by isotropic scattering (HG with r = 0). Fig. 3 compares our 
simulation results using this model with Antares measurements from June 
2000 [6] . The agreement verifies the validity of the model and our photon flux 
simulation for this case. 
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Figure 3. Residual time distributions at two distances from monochromatic flashes 
in water. The circles and stars are calibration data from ANTARES [6], while the 
solid curves are our simulation results using the measured water properties. The 
emission time profiles of the light sources (dashed curves) were measured in air, 
where scattering and absorption can be ignored. The distributions are normalised 
to unity at the peak value. 
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4-2 Modelling of ice and IceCube 's Antarctic glacial ice surveys 

The IceCube neutrino telescope, under construction deep in the glacial ice 
at the geographic South Pole, is planned to become a high-energy neutrino 
detector of 1 km 3 instrumented volume [14]. It is planned to have 80 strings, 
of which 22 have been deployed as of January 2007, each equipped with 60 
encapsulated photomultipliers evenly distributed over depths between 1450 m 
and 2450 m. IceCube builds on to the 19 strings of the Amanda array, which 
have been in operation since 2000. 

A detailed study of the properties of the deep South Pole glacial ice has been 
performed by the Amanda collaboration [8]. The glacial ice is very clear in 
the optical and near UV with absorption lengths of 20-120 m depending on 
wavelength. At wavelengths shorter than ~210 nm and longer than ~500 nm, 
absorption is dominated by the properties of pure ice, while in the intermedi- 
ate range absorption by impurities dominates. The effective scattering length 
is on average 25 m, less for shorter wavelengths. Both scattering and absorp- 
tion are strongly depth dependent and vary on all depth scales. The variations 
at depths exceeding 1450 m, where bubbles no longer exist, are explained by 
varying concentrations of insoluble dust deposits which correlate with changes 
in climatic conditions over geological time scales. By using physics motivated 
functional forms for the wavelength and depth dependences of scattering and 
absorption, the Amanda collaboration has elaborated a heterogeneous ice 
model [8] by investigating a large number of recorded light distributions gen- 
erated by in situ pulsed and steady light sources at different wavelengths. The 
resulting effective scattering and absorptions lengths, A e and A a as functions 
of wavelength, were averaged in 10 m depth intervals. 

Using the Amanda ice model parameters, we have generated simulated time 
distributions corresponding to two combinations of wavelength and light source- 
receiver positions, and compare them with experimental distributions in Fig. 4. 
The thick solid curves are our results when using the scattering and absorp- 
tion parameters fitted to these particular delay time distributions [8], with 
dashed curves representing two opposing deviations within the parameter un- 
certainty from these fits. The parameters fitted to the displayed distributions 
were A c = 27.6 m, A a = 20.5 m for the 532 nm curve, and A c = 22.6 m, 
A a = 82.0 m for the 470 nm curve, both with r = 0.94. The slight over- 
estimation in the tail in Fig. 4(b) is due to systematic uncertainties in the 
simulation of the LED emitter which lead to an imperfect description of the 
data in this particular case. The thin solid curves are simulation results with 
the heterogeneous ice model which is based on data from many source-receiver 
combinations at different depths and wavelengths. The differences between the 
two solid curves in each picture reflect the fact that the ice model parameters 
are averages over all fitted parameters in 10 m depth bins and the parameters 
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from individual fits are distributed around these averages. 
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Figure 4. Residual time distributions for two monochromatic pulsed light sources 
in deep glacial ice. In (a), light is recorded at a horizontal distance of 75 m from 
an isotropic laser source located at a depth of 1825 m. In (b), the detector is at 
a horizontal distance of 140 m from an upward-pointing LED emitter located at 
a depth of 1580 m. The Amanda calibration data are shown with Poissonian er- 
ror bars. The intrinsic timing widths of the light sources, less than 10 ns, were 
not included in the simulation. The thick solid curves are Photonics simulations 
using the scattering and absorption parameters fitted to these particular time dis- 
tributions, and the thin dashed curves represent two opposite parameter variations 
within the parameter uncertainties from the fits. The thin solid curves show our 
simulation results with the heterogeneous ice model [8]. 

5 Application to neutrino astronomy 

In neutrino astronomy, the universe is studied using high energy neutrinos as 
cosmic messengers. The neutrinos can be detected optically only after inter- 
acting with matter in the vicinity of the neutrino telescope and producing 
charged, Cherenkov light emitting particles like muons. Some of the emitted 
light is recorded by optical sensors distributed throughout the detector vol- 
ume. 

Ultra-relativistic muons are the primary channel through which high energy 
neutrinos are detected by optical neutrino telescopes. They are also the pri- 
mary background in the form of so-called atmospheric muons arising from high 
energy cosmic-ray interactions in the Earth's atmosphere. An extraterrestrial 
neutrino signal is distinguished from the background of neutrinos and muons 
created in the atmosphere mainly by differences in energy spectra and angu- 
lar distributions. It is therefore important to establish the particle energy and 
direction in every recorded event as accurately as possible. Our software con- 
tributes to this aim by providing means for detailed photon flux simulations, 
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using depth and wavelength dependent optical properties as established at a 
specific site. High-energy neutrino telescopes are typically recording data for 
several years while they are constructed by adding more optical sensors. Pho- 
ton flux tables generated by Photonics can cover arbitrarily large volumes 
and their use is therefore easily scalable to such growing sensor arrays. 

To further illustrate the utility of Photonics, we present in this section pho- 
ton propagation results for the inhomogeneous ice at the site of the IceCube 
neutrino telescope [8] . The flux in all figures is given per emitted photon, and is 
weighted with the angular and wavelength dependent acceptance of the optical 
detectors used by IceCube, so that the figures display the expected number of 
detections normalised to a 1 m 2 detection area in the direction of maximum 
optical module sensitivity. The photon flux tables were also convolved with a 
10 ns wide Gaussian to account for photomultiplier jitter. 

Fig. 5 shows the photon flux from a simulated infinitesimal electromagnetic 
cascade at 1730 m depth (Z s = in Amanda detector coordinates). Cascades 
are initiated in muon energy loss processes as well as in primary neutrino 
interactions. The light spectra of hadronic and electromagnetic cascades are 
Cherenkov in nature, but the light originates from many Cherenkov light emit- 
ting particles. The Cherenkov emission cone is slightly distorted [16] since not 
all emitting particles travel in parallel or at the same speed. The cascade in 
Fig. 5 is oriented toward the surface, at @ s = 135°, pointing at the upper 
left corner of the picture. A vertical slice through the photon flux containing 
the principal axis of the light source is displayed. The angular distribution 
of emitted photons is peaked at the Cherenkov angle, but the photon flux is 
smoothed out by scattering as it evolves through the ice. After 100 ns, we 
can still observe a peaked light distribution in the forward direction. At later 
times, the flux becomes more and more isotropic, to asymptotically resemble 
that of an isotropic flash. 

Since the IceCube sensors are pointing downward they detect upgoing light 
with a higher efficiency than downgoing light, which needs to be scattered to 
reach the photomultiplier photocathode. As a result, a point-like light source 
appears more upward. This can be seen in Fig. 5 where the direction of the 
light source appears to be oriented at an angle larger than 6 S = 135°. 

While the ice is very clear in the centre of the Amanda telescope, there are 
other depths where stronger scattering and absorption distort the photon flux. 
In the upper panel of Fig. 6 simulation results are shown for a setup analogous 
to that of Fig. 5 but with the shower origin 350 m deeper in the ice. This 
location is immediately below a region of strong scattering and absorption. The 
cascade direction is again S = 135°, but because of the particular medium 
properties in this region the event shape appears to be more isotropic and 
even resembles a downward pointing cascade. However, it does differ from 
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Figure 5. Light flux generated by a simulated electromagnetic cascade near the 
centre of the Amanda telescope. The upper panel shows a vertical slice of the photon 
flux (left) and the probability distribution (right) at t = 100 ns after light emission 
from a cascade at the origin and oriented toward the upper left at G s = 135°. The 
lower panel shows the same distributions at t = 300 ns. Note the different scales in 
the two panels. 



truly downward pointing events at the same location, shown in the lower panel 
of Fig. 6. In this particular case it is exceptionally hard to characterise the 
event correctly, but by using the correct heterogeneous medium description we 
improve the possibility to distinguish these cases. This is important for event 
simulation and reconstruction of parameters such as the zenith angle and the 
cascade energy. 

Inhomogeneities in the detector medium also strongly affect the optical topol- 
ogy of muon events. Fig. 7 shows the light distribution of a simulated muon 
moving upward through deep South Pole ice. At the front of the track, we 
observe a cross section of the unscattered Cherenkov wavefront, followed by 
a diffuse light cloud as the photons are scattered away from the geometri- 
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Figure 6. Light fluxes generated by two simulated electromagnetic cascades 350 m 
below the Amanda centre. This location is immediately below a region with higher 
dust concentration that causes stronger scattering and absorption. The flux is shown 
for both cascades at 100 ns and 300 ns after light emission. Because of the strong 
scattering and absorption just above the emission point, the flux from the upward 
pointing (0 S = 135°) source in the upper panel, (a) and (b), evolves from a dis- 
tribution peaked at the Cherenkov angle to a distribution similar to that from a 
downward pointing (0 S = 0°) source, shown in the lower panel, (c) and (d). 



cal Cherenkov cone. At depths with higher dust concentrations, photons are 
obstructed by scattering and absorbed before they can travel very far. This 
deforms the conical light front, which appears to be bent backwards. In the 
dusty region near z = —350 m, the photon flux is depleted and the surviving 
photons delayed by increased scattering. Muon track reconstruction is often 
strongly dependent on the earliest recorded photons, corresponding to the 
Cherenkov wavefront. Distortions in the wavefront, like in Fig. 7(b), can de- 
grade the reconstruction accuracy unless they are taken into account by the 
track fitting algorithms. 



22 



Fig. 8 shows a snapshot of the light field generated by a finite muon track with- 
out secondary interactions. The muon was created at the origin and propa- 
gated upward at 6 S = 135° until it decayed after 142 m. The probability 
density function for such a relatively short track approaches a shape similar 
to that of point-like cascades, making it hard to distinguish the two cases in 
an experimental situation with a limited number of points sampled by light 
sensors. PHOTONlcs-based simulations with realistic medium and light source 
descriptions allow experimentalists to isolate the differences in light profiles 
for these (and other) distinct cases and to develop appropriate reconstruction 
algorithms. 




horizontal position, x [m] 




horizontal position, x [m] 



(a) Ht) (b) /pdf(t) 

Figure 7. A snapshot of the light field produced by a muon which entered from 
below, at an angle O s = 135°. Inhomogeneities in the medium properties distort 
the smoothly arched light cone, as is most easily seen in the probability density 
function (b) in the particularly dusty region around z = —350 m which has stronger 
scattering and absorption. In (a), the flux is depleted in the dusty region. 
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Figure 8. Light distribution from a 142 m long muon track without secondary 
interactions. This muon was created at the origin and propagated upward at 
B s = 135° until it decayed at (x,z) = (—100,100). The figure shows a snapshot 
147 ns after the muon disappeared. Both the photon flux (a) and the probability 
density function (b) for such a comparably short track are similar to those produced 
by a point-like cascade. 
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6 Conclusion 



In this paper we have presented the concepts and methods which combine into 
the Photonics software package. We have explained how the program can 
be used for calculating and tabulating light distributions around a stationary 
or moving source, as a function of time and space in scattering and absorbing 
heterogeneous media. The light distributions obtained from our Monte Carlo 
simulation agree well with observations of calibration light sources in deep 
sea water and glacial ice surveys. In the last section it is demonstrated how 
Photonics can be used to model how optical inhomogeneities of the Antarctic 
ice at the location of the IceCube neutrino telescope distort the footprints of 
elementary particle interactions. 
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